############## 
############## 
############## Models
remove(list = ls())

base::library(conflicted)
base::library(tidyverse)
conflict_prefer("filter","dplyr")
base::library(ggplot2)
base::library(dplyr)
base::library(here)
conflict_prefer("here", "here")
base::library(systemfit)
base::library(stargazer)
base::library(DataCombine)
base::library(fixest)
base::library(sandwich)
base::library(stargazer)
base::library(ggeffects)

load(here("Data", "Model_Data.RData"))
glimpse(final_data)

summary(final_data$death_rate)
table(final_data$year)

final_data <- final_data %>% 
  mutate(time_trend = ifelse(year == 2004, 1, NA),
         time_trend = ifelse(year == 2008, 2, time_trend),
         time_trend = ifelse(year == 2012, 3, time_trend),
         time_trend = ifelse(year == 2016, 4, time_trend),
         time_trend = ifelse(year == 2020, 5, time_trend)) 

final_data <- final_data %>% 
  mutate(fouryear_chagen = death_rate-lastelec_death_rate) %>% 
  mutate(pc_4change_death = 100*(fouryear_chagen)/lastelec_death_rate)

glimpse(final_data)

eq1_d <- dem_rep ~ pc_4change_death + unemp_rate + 
  log(defl_pcincome) + 
  share_black + share_hisp + share_asian + share_young + share_elder + 
  log1p(share_college) + log1p(share_manuf) +
  log(tot_pop) + as.factor(rural_urban) + state_partiesdiff +
  state + rep_inc + incumbent_cand + lag_dem_rep + time_trend
eq2_d <- abs_rep ~ pc_4change_death + unemp_rate + 
  log(defl_pcincome) + 
  share_black + share_hisp + share_asian + share_young + share_elder + 
  log1p(share_college) + log1p(share_manuf) +
  log(tot_pop) + as.factor(rural_urban) + state_partiesdiff +
  state + rep_inc + incumbent_cand + lag_abs_rep + time_trend
Syst_death <- list(eq1_d, eq2_d)
model_d <- systemfit(Syst_death, method = "SUR", data = final_data)
summary(model_d)



ols_model1 <- lm(rep_share2 ~ pc_4change_death + unemp_rate + log(defl_pcincome) +
                   share_black + share_hisp + share_asian + share_young + share_elder + 
                   log1p(share_college) + log1p(share_manuf) +
                   log(tot_pop) + as.factor(rural_urban) + state + state_partiesdiff + 
                   rep_inc + incumbent_cand + time_trend + lag_rep_share2,
                 final_data); summary(ols_model1)


fe_model <- feols(rep_share2 ~ pc_4change_death + unemp_rate + log(defl_pcincome) +
                    share_black + share_hisp + share_asian + share_young + share_elder + 
                    log1p(share_college) + log1p(share_manuf) +
                    log(tot_pop) + as.factor(rural_urban) + state + state_partiesdiff + 
                    rep_inc + incumbent_cand + time_trend + lag_rep_share2 | county,
                  data = final_data); summary(fe_model)


summary(ols_model1)

ols_model2 <- coeftest(ols_model1, vcov = vcovCL, cluster = ~ county)
ols_model2


stargazer(ols_model1, ols_model2, ols_model1,
          type = "latex",
          title = "OLS and Fixed Effect Estimates",
          dep.var.labels = c("Republican Vote Share", "Republican Vote Share", "Republican Vote Share"), 
          column.labels = c("State Fixed Effecs", 
                            "Clustered SE by County", 
                            "County Fixed Effects"),
          no.space = FALSE, single.row = T,
          covariate.labels = c("Perc. Change in Overdose Death Rate",
                               "Local Unemployment Rate",
                               "Log PC Income",
                               "Share Black",
                               "Share Hispanic",
                               "Share Asian",
                               "Share Young Voters",
                               "Share Elder Voters",
                               "Log Share College Degree",
                               "Log Share Manufacture",
                               "Log Population",
                               "Close Election (State)",
                               "Republican Presidency",
                               "Incumbent Candidate",
                               "Lagged Republican Vote Share",
                               "Constant"),
          omit = c("rural_urban", "state"),
          omit.labels = c("Rural-Urban Code",
                          "State Fixed Effects"))




summary(final_data$pc_4change_death)
ols_model2

model_rep2 <- fe_model %>% 
  ggpredict(terms = c("pc_4change_death[18, 20, 22, 24, 26, 28, 30, 32, 34, 36]"), 
           ci.lvl = 0.90)

glimpse(model_rep2)

model_drugs2 <- data.frame(model_rep2)
glimpse(model_drugs2)


change_pic <- model_drugs2 %>% 
  ggplot(aes()) + 
  geom_pointrange(aes(x = x, y = predicted, 
                      ymin = conf.low, ymax = conf.high),
                  lwd = .85, 
                  color = "darkslategray"
  ) + 
  theme_light() + 
  labs(x = "% Change in Death Rates", y = "Predicted Republican Vote Share",
       title = NULL) +
  ggtitle(NULL) +
  scale_x_continuous(breaks = c(20, 25, 30, 35), lim = c(18, 36)) +
  scale_y_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1), lim = c(0, 1)) +
  theme(axis.text.x = element_text(hjust = 0.5, size = 11),
        axis.text.y = element_text(size = 11),
        axis.title.y = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        legend.position = "bottom",
        legend.direction = "horizontal",
        plot.title = element_text(size = 14))

change_pic








